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Abstract 


We  represent  the  standard  ramp  filter  operator  of  the  filtered  back-projection  (FBP)  recon¬ 
struction  in  different  bases  composed  of  Haar  and  Daubechies  compactly  supported  wavelets. 

The  resulting  multiscale  representation  of  the  ramp  filter  matrix  operator  is  approximately  diag¬ 
onal.  The  accuracy  of  this  diagonal  approximation  becomes  better  as  wavelets  with  larger  num¬ 
ber  of  vanishing  moments  are  used.  This  wavelet-based  representation  enables  us  to  formulate  a 
multiscale  tomographic  reconstruction  technique  wherein  the  object  is  reconstructed  at  multiple 
scales  or  resolutions.  A  complete  reconstruction  is  obtained  by  combining  the  reconstructions  at 
different  scales.  Our  multiscale  reconstruction  technique  has  the  same  computational  complexity 
as  the  FBP  reconstruction  method.  It  differs  from  other  multiscale  reconstruction  techniques  in 
that  1)  the  object  is  defined  through  a  multiscale  transformation  of  the  projection  domain,  and 
2)  we  explicitly  account  for  noise  in  the  projection  data  by  calculating  maximum  aposteriori 
probability  (MAP)  multiscale  reconstruction  estimates  based  on  a  chosen  fractal  prior  on  the 
multiscale  object  coefficients.  The  computational  complexity  of  this  MAP  solution  is  also  the 
same  as  that  of  the  FBP  reconstruction.  This  is  in  contrast  to  commonly  used  methods  of  sta¬ 
tistical  regularization  which  result  in  computationally  intensive  optimization  algorithms.  The 
framework  for  multiscale  reconstruction  presented  here  can  find  application  in  object  feature 
recognition  directly  from  projection  data,  and  regularization  of  imaging  problems  where  the 
projection  data  are  noisy. 

Key  words;  multiresoluiion  reconstruction,  wavelets,  tomography,  stochastic  models. 
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1  Introduction 


In  this  work  we  present  a  multiresolution  approach  to  the  problem  of  reconstructing  an  image  from 
tomographic  projections.  In  general,  a  multiresolution  framework  for  tomographic  reconstruction 
may  be  natural  or  desirable  for  a  variety  of  reasons.  First,  the  data  under  consideration  may 
be  naturally  acquired  at  multiple  resolutions,  e.g.  if  data  from  detectors  of  differing  resolutions 
are  used.  In  addition,  the  phenomenon  may  itself  be  naturally  multiscale.  For  example,  in  the 
medical  field  self-similcir  or  fractal  models  have  been  effectively  used  for  the  liver  and  lung  [10-12]. 
Furthermore,  it  may  be  that,  even  if  the  data  and  phenomenon  are  not  multiscale,  our  ultimate 
objectives  are  multiresolution  in  some  way.  For  example,  even  though  our  data  may  be  acquired  at 
a  fine  level  we  may  actually  only  care  about  aggregate  or  coarse  scale  quantities  of  the  field.  Such 
is  often  the  case  in  ocean  acoustic  tomography  or  frmctional  medical  imaging.  Conversely,  if  we  are 
only  interested  in  imaging  high  frequency  details  within  the  object  (for  example,  boundaries),  then 
we  could  directly  obtain  these  features  by  extracting  only  the  finer  scade  information  in  the  data.  Or, 
indeed,  it  may  be  that  we  want  to  use  different  resolutions  in  different  areas  -  e.g.  in  nondestructive 
evaluation  of  aircraft  we  may  want  to  look  for  general  corrosion  over  an  entire  plane,  but  focus 
attention  on  certain  suspect  rivets  to  look  for  cracks.  Using  conventional  techniques  we  would 
first  have  to  reconstruct  the  entire  field  and  then  use  post-processing  to  extract  such  features.  A 
final  compelling  motivation  for  the  use  of  multiresolution  methods  in  estimation  emd  reconstruction 
problems,  is  that  they  generally  lead  to  extremely  efficient  algorithms,  as  in  [13]. 

The  conventional,  and  most  commonly  used,  method  for  image  reconstruction  from  tomographic 
projections  is  the  Filtered  Back-Projection  (FBP)  reconstruction  technique  [1].  In  the  standard 
FBP  reconstruction  as  applied  to  a  complete  set  of  noiseless  projections^  the  projection  data  at  each 
cingle  are  first  filtered  by  a  high-pass  “ramp”  filter  and  then  back-projected.  In  this  paper,  we  work 
in  a  different,  multiscale  transform  space.  The  matrix  representation  of  the  resulting  multiscale 
filtering  operator  is  approximately  diagonal.  This  enables  us  to  formulate  an  efficient  multiscale 
tomographic  reconstruction  technique  that  has  the  same  computational  complexity  as  that  of  the 
FBP  reconstruction  method.  Perhaps  more  significamtly,  however,  the  different  scale  components 
of  our  proposed  multiscale  reconstruction  method  induce  a  corresponding  mTiltisczde  representation 
of  the  underlying  object  and,  in  particulair,  provide  estimates  of  (and  thus  information  about)  the 
field  or  object  at  a  variety  of  resolutions  at  no  additional  cost.  This  provides  a  natural  framework 
for  explicitly  assessing  the  resolution- accuracy  tradeoff  which  is  critical  in  the  case  of  noisy  or 
incomplete  data. 

Noisy  imaging  problems  arise  in  a  variety  of  contexts  (e.g.  low  dose  medical  imaging,  oceanog¬ 
raphy,  and  in  several  applications  of  nondestructive  testing  of  materials)  and  in  such  cases  standard 
techniques  such  as  FBP  often  yield  unacceptable  results.  These  situations  generally  reflect  the  fact 
that  more  degrees  of  freedom  are  being  sought  than  are  really  supported  by  the  data  and  hence 
some  form  of  regularization  is  required.  Conventionally,  this  problem  of  reconstruction  from  noisy 
projection  data  is  regularized  by  one  of  the  following  two  techniques.  First,  the  FBP  ramp  filter 
may  be  rolled  off  at  high  frequencies  thus  attenuating  high  frequency  noise  at  the  expense  of  not 
reconstructing  the  fine  sczde  features  in  the  object  [7,8].  This  results  in  a  fast,  though  ad  hoc, 
method  for  regularization.  The  other  common  method  for  regularization  is  to  solve  for  a  maximum 
aposteriori  probability  (MAP)  estimate  of  the  object  based  on  a  2-D  (spatial)  Markov  random  field 
(MRF)  prior  model  [25,26].  This  results  in  a  statistically  regularized  reconstruction  which  allows 

‘According  to  Llacer  [3],  “a  complete  data  set  could  be  described  as  sufficient  number  of  line  projections  at  a 
sufficient  number  of  angular  increments  such  that  enough  independent  measurements  are  made  to  allow  the  image 
reconstruction  of  a  complete  bounded  region.” 
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the  inclusion  of  prior  knowledge  in  a  systematic  way,  but  leads  to  optimization  problems  which  are 
extremely  computationadly  intensive.  In  contrast  to  these  methods,  we  are  able  to  extend  our  mul¬ 
tiscale  reconstruction  technique  in  the  case  of  noisy  projections  to  obtain  a  multiscale  MAP  object 
estimate  which,  while  retaining  all  of  the  advantages  of  statistically-based  approaches,  is  obtained 
with  the  same  computational  complexity  as  the  FBP  reconstruction.  We  do  this  by  realizing  that 
for  ill-posed  problems  the  lower  resolution  (i.e.  the  coarser  scale)  reconstructions  are  often  more 
reliable  than  their  higher  resolution  comterparts  and  by  capturing  such  intuition  in  prior  statistical 
models  constructed  directly  in  the  multiscale  domain.  Similar  to  the  noise-free  case,  we  also  obtain 
these  MAP  estimates  at  multiple  scales,  essentially  for  free. 

Finally,  the  FBP  reconstruction  method  is  not  suitable  for  imaging  problems  where  the  projec¬ 
tion  data  are  incomplete  (limited  angle  and/or  truncated  projections)  [2,3].  These  problems  are 
encountered  in  metny  applications  in  medicine,  non-destructive  testing,  oceanography  and  surveil- 
lemce.  Though  the  work  presented  here  focuses  on  the  case  of  complete  data,  our  wavelet-based 
multiscale  framework  has  the  potential  of  overcoming  this  limitation  and  we  briefly  discuss  such 
possibilities  in  the  conclusion  to  the  paper.  We  also  refer  the  reader  to  a  subsequent  paper  [20] 
where  we  consider  an  extension  to  the  incomplete  data  case  based  on  a  similar  multiscale  framework. 

Wavelets  have  been  recently  applied  to  tomography  by  other  researchers  as  well.  Peyrin  et 
al  [15]  have  shown  that  the  back-projection  of  ramp-flltered  and  wavelet-transformed  projection 
data  corresponds  to  a  2-D  wavelet  decomposition  of  the  original  object.  Their  method  differs  from 
ours  in  several  ways.  First,  the  work  in  [15]  does  not  deal  with  noise  in  the  projection  data.  In 
contrast,  our  frcunework  allows  for  the  solution  of  statistically  regularized  problems  at  no  additional 
cost  when  the  projection  data  are  noisy.  Second,  in  [15],  the  object  is  represented  in  the  original 
spatial  domain  by  a  2-D  wavelet  basis,  the  expansion  coefficients  of  which  are  then  obtained  from 
the  projection  data.  Instead  of  this,  we  represent  the  object  in  the  projection  domain  by  expemding 
the  FBP  basis  functions  (i.e.  strips)  in  a  1-D  wavelet  basis.  This  has  the  advantage  that  our 
multiscale  basis  representation  of  the  object  is  closer  to  the  measurement  domain  than  the  multiscale 
representation  in  [15].  One  consequence  is  that  our  algorithms  for  multiscale  reconstruction  are  no 
more  complex  tham  the  FBP  method.  Another,  is  that  our  framework  also  allows  for  the  simple 
zmd  efficient  solution  of  statistically  regularized  problems  at  no  additional  cost  when  the  projection 
data  are  noisy.  Scdiiner  zind  Yagle  use  the  wavelet  transform  to  perform  spatially-varying  Altering 
by  reducing  the  noise  energy  in  the  reconstructed  image  over  regions  where  high-resolution  features 
are  not  present  [16].  They  also  apply  wavelet  based  reconstruction  to  the  limited  eingle  tomography 
problem  by  assuming  approximate  a  priori  knowledge  about  the  edges  in  the  object  that  lie  parallel 
to  the  missing  views  [17].  Again  as  in  [15],  in  [16,17]  the  object  is  represented  in  the  original  spatial 
domeiin  by  a  2-D  wavelet  basis  which  is  much  different  than  our  multiscale  object  representation. 
DeStefano  and  Olson  [18],  and  Berenstein  and  Walnut  [19]  have  also  used  wavelets  for  tomographic 
reconstruction  problems,  in  particular  to  localize  the  Radon  transform  in  even  dimensions.  Through 
this  localization  the  radiation  exposure  can  be  reduced  when  a  local  region  of  the  object  is  to  be 
imaged.  The  work  in  [18]  and  [19]  does  not  provide  a  framework  for  multiscale  reconstruction, 
however,  which  is  the  central  theme  here.  In  addition,  our  reconstruction  procedure  also  locjilizes 
the  Radon  transform,  though  we  do  not  stress  this  particular  application  in  this  work. 

The  paper  is  organized  as  follows.  Section  2  contains  preliminaries.  In  Section  2.1  we  de¬ 
scribe  the  standard  tomographic  reconstruction  problem  and  in  Section  2.2  we  describe  the  FBP 
reconstruction  technique.  We  outline  the  theory  of  1-D  multiscale  decomposition  in  Section  2.3. 
In  Section  3,  we  develop  the  theory  behind  our  wavelet-based  multiscale  reconstruction  method 
starting  from  the  FBP  object  representation.  In  Section  4  we  btuld  oh  this  framework  to  provide  a 
fast  method  for  obtaining  MAP  regularized  reconstructions  from  noisy  data.  The  conclusions  are 
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:  projection  at  angle  1  (k=l) 


Figure  1;  The  projection  measurements  for  an  object,  /  (shaded),  at  two  different  angular  positions 
(fc  =  1  and  k  =  2  respectively).  The  number  of  parallel  strip  integrals  in  each  angular  projection, 
N,,  is  8  in  this  case.  Three  basis  functions,  Tii,ri8,T28,  which  are  the  indicator  functions  of  the 
corresponding  strips,  are  also  shown. 


presented  in  Section  5.  Appendix  A  summarizes  the  mathematical  notation  used  throughout  this 
paper.  Appendix  B  contains  cert2un  technical  details. 

2  Preliminaries 

2.1  The  Tomographic  Reconstruction  Problem 

In  tomography  the  goal  is  to  reconstruct  an  object  or  a  field,  /,  from  Une-integral  projection  data  [1]. 
For  a  parallel-beam  imaging  geometry,  the  projection  data  consists  of  parallel,  non-overlapping  strip 
integrals  through  the  object  at  veu'ious  angles  (refer  to  Figure  1).  Each  emgular  position  corresponds 
to  a  specific  source- detector  orientation.  Suppose  we  have  Ng  uniformly  spaced  angular  positions 
between  0°  and  180®  and  N,  parallel  strip  integrals  at  each  angular  position.  Let  us  label  the 
observation  corresponding  to  projection  I  at  angular  position  k  by  yk{^),  where  k  =  l,...,Ng  emd 
I  =  1,. .  .,Ng.  Furthermore,  let  Tki  be  the  indicator  function  of  the  strip  integral  corresponding  to 
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this  observation  so  that  Tju  has  value  one  within  that  strip  and  zero  otherwise.  Given  this  notation, 
yk{i)=JJ  f{u,v)Tktiu,v)dudv  k  =  £  =  1,. .  .,Nt  (1) 

n 

where  (u,  v)  are  the  usual  rectangular  spatial  coordinates  and  the  integration  is  carried  over  a  region 
of  interest 

Due  to  practical  considerations,  we  have  to  work  with  a  discretized  version  of  (1).  By  using 
standard  discretization  techniques  (see  for  example  [21]),  the  projection  data  at  angle  k  is  given  by 

Vk  =  Tkf  (2) 

where  Tk  is  Em  iV,  x  N,  matrix  representing  {Tki{u,  v)-,  £  =  1, iV,}  and  /  is  em  x  1  vector 
representing  f{u,v)  on  an  i\r,  x  Ng  square  pixel  lattice,  and  yk  is  the  corresponding  vector  of 
measurements  yk{£)-  Thus  row  £  of  Tk  is  the  (discrete)  representation  of  the  strip  function  Tki{u,  t») 
and  the  inner  product  of  /  with  this  strip  yields  the  data  contained  in  the  corresponding  entry  of  y*. 
The  tomographic  reconstruction  problem  then  reduces  to  finding  an  estimate  /  of  the  discretized 
object  f  given  the  projection  data  contained  in  the  {yk’,  k  =  1,. . N$}. 

2.2  The  Filtered  Back-Projection  Reconstruction  Technique 

The  filtered  back-projection  (FBP)  reconstruction  technique  is  the  most  commonly  used  method  for 
image  reconstruction  from  tomographic  data.  It  is  based  directly  on  the  Radon  inversion  formula 
which  is  valid  (i.e.  yields  exact  reconstructions)  only  when  a  continuum  of  noise-free  line  integreil 
projections  from  aU  angles  are  used  [Ij.  In  practice,  as  indicated  in  (1),  we  only  have  access  to 
sampled  projection  data  which  Eire  collected  using  strips  of  finite  width.  In  this  case,  the  qucdity 
of  the  FBP  reconstruction  is  a  function  of  the  quality  and  fineness  of  the  corresponding  projection 
data  used.  We  refer  the  reader  to  [23,24]  for  details  on  sampling  requirements  for  the  Radon 
trsmsform.  In  this  work  we  assume  that  we  sample  finely  enough  to  produce  good  reconstruction 
in  the  noiseless  case. 

In  the  FBP  reconstruction,  the  object  is  expanded  in  a  non-orthogonal  basis  that  is  closely 
related  to  the  data  acquisition  process  Eind  the  coef5.cients  of  this  fexpEuision  are  then  found  from 
locEd  processing  of  the  data  in  each  angular  projection.  In  paxticulEir,  the  estimated  object  is 
represented  as  a  linear  combination  of  the  SEime  functions  Tki{u,  v)  along  which  the  projection  data 
are  collected.  Similar  to  (2),  a  discretized  version  of  this  representation  may  be  obtEuned  as: 

k=l 

where  the  N,  vector  i*.  contains  the  object  coefficient  set  at  Euigle  k.  Note  that  (3)  can  be  interpreted 
as  the  back-projection  operation  [1]  where  the  object  coefficients  Xk  are  back-projected  along  the 
basis  fimctions  Tk  at  each  angle  k  and  then  the  contributions  from  aU  Ng  Eingles  are  added  to  get 
the  overEiU  reconstruction  /. 

To  complete  the  reconstruction  the  coefficients  Xk  must  now  be  determined.  The  standard  FBP 
method  calculates  them  for  each  angle  k  according  to  the  Radon  inversion  formula  by  filtering  the 
projection  data  yk  at  that  pEirticular  angle  with  a  ramp  filter  [1].  Thus,  for  a  fixed  angle  k: 

Xk  =  Ryk  (4) 
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where  the  matrix  R  captures  this  ramp-filtering  operation.  Thus  (3)  and  (4)  together  represent  the 
two  operations  used  in  the  standard  FBP  reconstruction. 

In  Section  3  we  apply  a  1-D  multiscale  change  of  basis  to  the  coefficients  Xk  and  observations  yu 
using  the  wavelet  transform  [4,22].  One  effect  of  this  multiscale  tramsformation  will  be  to  “compress” 
(sparsify)  the  filtering  operator.  Beyond  simple  compression  of  this  operator,  however,  our  method 
results  in  an  associated  multiscade  transformation  of  the  basis  fimctions  contained  in  T^,  aind  thus 
leads  naturally  to  a  framework  for  the  reconstruction  of  objects  at  multiple  resolutions,  and  hence 
the  extraction  of  mformation  at  multiple  resolutions.  A  key  point  in  our  multiscale  reconstruction 
method  is  that,  as  opposed  to  what  is  done  in  other  multiscale  reconstruction  techniques  (for 
example,  [15]),  we  do  not  directly  expand  the  object  estimate  (i.e.  /)  in  a  spatial  2-D  wavelet  basis 
but  rather  we  expand  the  FBP  basis  functions  {2*,}  in  a  1-D  wavelet  basis  which  then  induces  a 
corresponding  2-D  multiscale  object  representation.  The  multiscale  versions  of  the  object  expansion 
coefficients,  {**},  “live”  in  the  strip  integral  (i.e.  the  projection)  domain  rather  than  in  the  original 
object  or  spatial  domain.  Thus,  as  we  have  said,  our  multiscale  basis  representation  of  the  object  is 
closer  to  the  measurement  domain  than  the  multiscale  representations  in  previous  work,  resulting 
in  very  efficient  algorithms. 

2.3  1-D  Wavelet  Transform  Based  Multiscale  Decomposition 

Here  we  present  a  brief  summary  of  the  wavelet-based  multiscale  decomposition  of  1-D  functions. 
We  do  not  intend  this  as  a  complete  treatment  of  the  topic  and  intentionally  suppress  many  details. 
The  interested  reader  is  referred  to  any  of  the  many  papers  devoted  to  this  topic,  e.g.  the  excellent 
paper  [22].  A  multiresolution  dyadic  decomposition  of  a  discrete  1-dimensional  signal  x(n)  of  length 
2^  is  a  series  of  approximations  ®("*^(n)  of  that  signal  at  finer  and  finer  resolutions  (increasing  m) 
with  dyadicaUy  increasing  complexity  or  length  and  with  the  approximation  at  the  finest  scale 
equaling  the  signal  itself,  i.e.  x^^\n)  =  x{n).  By  considering  the  incremental  detail  added  in 
going  from  the  m-th  scale  approximation  to  the  (m  -f  l)-st  we  arrive  at  the  wavelet  transform.  In 
particular,  if  is  the  vector  contcdning  the  sequence  E^”*)(n)  and  is  the  corresponding  vector 
of  detail  added  in  proceeding  to  the  next  finer  scale,  then  one  can  show  [27]  that  the  evolution  of 
the  approximations  in  scale  satisfies  an  equation  of  the  form: 

^  J.(m)  ^(m)  (5) 

where  L{m)  and  H (m)  are  matrices  (linear  transformations)  which  depend  on  the  particular  wavelet 
chosen  and  are  fax  from  arbitrary  and  L^{m)  and  H^{m)  denote  their  transposes  (i.e.  their  ad- 
joints).  The  operators  L^{m)  and  H^{m)  serve  to  interpolate  the  “low”  and  “high”  frequency  (i.e. 
approximation  and  detail)  information,  respectively,  at  one  scale  up  to  the  next  finer  scale.  The 
2’”-vector  containing  the  information  added  in  going  from  scale  m  to  (m  -i  1),  is  composed  of 
the  wavelet  transform  coefficients  at  scale  m  and  (5)  is  termed  the  wavelet  synthesis  equation. 

Starting  from  a  “coarsest”  approximation  (usually  t^lken  to  be  the  average  value  of  the 
signal)  then,  it  is  possible  to  recursively  and  efficiently  construct  the  different  scale  approximations 
through  (5)  by  using  the  complete  set  of  wavelet  coefficient  vectors  This  layered  construction 

is  shown  graphically  in  Figure  2a,  where  our  approximations  are  refined  though  the  addition  of 
finer  and  finer  levels  of  detail  as  we  go  from  right  to  left  until  the  desired  scale  of  approximation 
is  achieved.  In  particular,  the  original  signal  x  is  obtained  by  adding  all  the  interscale  detail 
components  to  the  initial  approximation  For  a  given  signal  x  the  complete  set  of  these 
elements  uniquely  captures  the  signal  and  thus  corresponds  to  a  simple  change  of  basis.  In  addition, 
note  from  Figure  2  that  the  intermediate  approximation  at  scale  m  is  generated  using  only 
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Figure  2:  (a)  Tree  diagram  for  wavelet  transform  synthesis.  We  start  from  a  coarsest  approximation 
on  the  right  and  progressively  add  finer  levels  of  detail  as  we  proceed  to  the  left,  thus 
refining  the  original  approximation  to  the  signal.  The  original  (finest  scale)  sequence  is  obtained 
as  the  final  output  on  the  left,  (b)  Tree  diagram  for  wavelet  transform  anedysis.  Starting  from  a 
finest  level  signal  x  in  the  left  we  recursively  peel  off  layers  of  detail  as  we  proceed  to  the  right 

and  the  next  coauser  scale  representation 

the  corresponding  subset  of  the  complete  wavelet  coefficient  set  (e.g.  to  obtain  we  use  oidy 
a;(°),  and  The  representation  of  this  intermediate  approximation  at  the  original  finest 

scade  can  be  found  by  repeated  interpolation  of  the  information  in  through  the  application 
of  m'  >  m.  This  interpolation  up  to  the  finest  scade  corresponds  to  effectively  assuming 

that  additional,  finer  scale,  detadl  components  m!  >m  are  zero  in  our  representation  of  the 

signal.  It  is  such  intermediate  scale  approximations  aind  the  detail  necessary  to  go  between  them 
that  give  the  wavelet  trainsform  its  natural  multiscade  interpretation,  and  indeed  we  exploit  such 
interpretations  in  Sections  3  and  4  to  obtain  induced  multiscale  ofiject  representations. 

Beyond  the  recursive  computation  of  the  approximations,  it  is  also  possible  to  compute  the 
components  of  the  decomposition  itself  (i.e.  the  wavelet  coefficients)  recursively  by  exploiting  the 
same  multiscale  structure.  In  particular,  as  shown  in  [27]  the  wavelet  coefficient  vectors  (and 
®(°))  can  be  obtained  from  the  following  recursion  defining  the  wavelet  analysis  equations,  which  is 
illustrated  in  Figure  2b: 

=  L{m)  =  H{m)  (6) 

where  L{m)  and  H{m)  are  the  same  operators  defined  in  connection  with  (5).  The  operators 
i(m)  md  Hijn)  correspond  roughly  to  low  and  high  pass  filters  followed  by  downsampling  by  a 
factor  of  2,  respectively.  The  figure  shows  how  these  wavelet  coefficient  vectors  at  each  scale  are 
obtained  by  “peeling  off”  successive  layers  of  detail  as  we  proceed  from  finer  to  coarser  scales  (left 
to  right  in  the  figure).  This  recursive  structure  yields  algorithms  for  computation  of  the  wavelet 
transform  coefficients  that  are  extremely  efficient.  For  convenience  in  the  development  to  follow,  we 
win  capture  the  overall  operation  which  takes  a  vector  x  containing  a  discrete  signal  to  the  vector 
£  containing  all  of  its  corresponding  wavelet  transform  elements  and  by  the  matrix 
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operator  W  as  follows: 


Since  the  transform  is  invertible  and  the  wavelet  basis  functions  zire  orthonormal,  it  follows  that 
W~^  exists  and  further  that  W  is  a  unitary  matrix,  i.e.  that  W~^  =  W^.  Prom  the  above  discussion, 
the  matrix  W  captures  the  operation  of  the  operators  L(m)  and  Him),  and  thus  depends  on  the 
underlying  chosen  wavelet.  In  our  work  in  this  paper,  in  addition  to  the  Haar  wavelet  we  wiU  use 
wavelets  from  an  especially  popular  family  of  these  functions  due  to  Daubechies  [4],  the  sepsuate 
elements  of  which  are  denoted  where  n  is  an  indication  of  the  support  size  of  the  corresponding 
filters  contained  in  L{m)  and  H{m).  Finally,  since  our  signals  are  of  finite  length,  we  need  to  deal 
with  the  edge  effects  which  occur  at  the  ends  of  the  interval  in  the  wavelet  transform.  While  there 
are  a  variety  of  ways  in  which  to  do  this,  such  as  modifying  the  wavelet  functions  at  the  ends  of  the 
interval  in  order  to  provide  an  orthogonal  decomposition  over  the  intervad  [28],  we  have  chosen  here 
to  use  one  of  the  most  commonly  used  methods,  namely  that  of  cyclically  wrapping  the  interval 
which  induces  a  circulant  structure  in  L{m)  and  H{m)  [5,22].  While  this  does  introduce  some  edge 
effects,  these  are  of  negligible  importance  for  the  objectives  and  issues  we  wish  to  emphasize  and 
explore  and  for  the  appUcations  considered  here.  Further,  the  methods  we  describe  can  be  readily 
adapted  to  other  approaches  for  dealing  with  edge  effects  as  in  [28]  and  the  references  contained 
therein. 

As  noted  above,  the  intermediate  approximations  and  their  finest  scale  representation 
may  be  obtained  by  using  only  part  of  the  full  wavelet  coefficient  set  during  synthesis,  effectively 
assuming  the  finer  scale  detail  components  are  zero.  For  convenience  in  the  discussion  to  follow 
we  capture  this  partial  zeroing  operation  in  the  matrix  operator  Aim),  which  nulls  the  upper 
N  —  m  subvectors  of  the  overall  wavelet  vector  ^  and  thus  retaiins  only  the  information  necessary 
to  construct  the  approximation  x^*")  at  scale  m: 

Aim)  =  block  diag  [0(2iv^_2'»),  -^(2™)]  (8) 

where  Op  is  a  p  X  p  matrix  of  zeros  and  Iq  is  a  q  X  q  identity  matrix.  Also  it  wifi,  prove  convenient 
to  define  a  similar  matrix  operator  Dim),  that  retains  only  the  information  in  ^  pertaining  to  the 
detail  component  at  scale  m  by  zeroing  all  but  the  sub- vector  corresponding  to 

Dim)  =  block  diag  j^0(2w^-2"‘+i)>  -^(2™))  0(2™)]  (9) 

Finally,  with  these  definitions  note  that  we  have  the  following  scale  recmrsive  relationship  for  the 
partially  zeroed  vectors,  in  the  spirit  of  (5): 

A(m+i)  ^  =  Aim)  ^  -f  Dim)  ^  (10) 

3  The  Multiscale  Reconstruction  Technique 

In  this  section  we  derive  our  1-D  wavelet-based  multiscale  reconstruction  technique.  We  stMt  by 
applying  a  wavelet-derived  multiscale  change  of  basis  W  to  the  FBP  object  coefficients  x*,  which 
will  induce  a  natural  multiresolution  object  representation.  We  then  show  how  the  coefficients  of  om 
new  multiscsde  representation  can  be  computed  directly  from  corresponding  multiscale  versions  of 
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the  data,  in  the  same  way  that  Xk  is  computed  directly  from  yk  in  the  standard  FBP  method.  Taken 
together  these  two  components  define  a  multiscale  reconstruction  zilgorithm,  analogous  in  structure 
to  the  FBP  method.  An  important  point  is  that  our  approach  does  not  start  with  a  decomposition 
of  the  object  in  a  2-D  wavelet  basis  and  attempt  to  then  find  the  resulting  coefficients,  but  rather 
works  directly  in  the  projection  domain.  The  multiscale  nature  of  our  object  representation  in 
the  2-D  or  spatial  domain  arises  naturally  from  the  original  FBP  definitions  and  our  multiscale 
decomposition  of  the  art,  and  thus  we  retain  the  simplicity  and  efficiency  of  this  poptdar  method. 

Multiscale  Object  Representation 

We  start  by  applying  a  multiscale  change  of  basis,  as  defined  by  the  matrix  W  in  Section  2.3,  to 
the  original  set  of  object  coefficients  at  each  angle  k  to  obtain  an  equivalent  set  of  multiscale 
object  coefficients  as  follows: 

a  =  Wxk  (11) 

Thus,  for  a  given  choice  of  wavelet  defining  W,  the  vector  contains  the  corresponding  wavelet 
coefficients  and  coarsest  level  approximation  (i.e.  the  average)  associated  with  Xk  and  thus  forms  a 
multiresolution  representation  of  this  signal.  More  importantly,  by  reflecting  this  change  of  basis 
into  the  original  FBP  object  representation  (3),  we  naturally  induce  a  corresponding  multiscale 
representation  of  the  object  through  the  creation  of  a  corresponding  set  of  transformed  multiscale 
basis  functions.  In  particular,  substituting  (11)  into  (3)  we  obtain: 

Ng 

/  =  i:  [Wx,)  a  ■£  (12) 

*=1  Jfe=l 

where  Tk  =  W  Tk,  is  now  a  matrix  representing  the  transformed,  multiscale  basis  functions  at  angle 
k. 

Before  proceeding,  let  us  consider  these  transformed  bases  functions  contained  in  Tk  in  more 
detail.  Recall  from  Section  2.1  that  the  rows  of  Tk  are  composed  of  the  (discretized)  original  strip 
basis  functions  at  angle  k  along  which  the  data  were  collected,  c.f.  (2).  Similarly  the  rows  of  the 
transformed  matrix  Tk  will  contain  the  corresponding  (discretized)  multiscale  object  basis  functions 
at  angle  k.  The  wavelet  transform  operator  matrix  W,  acting  identically  on  each  column  of  Tk, 
will  thus  form  the  new  multiscale  basis  functions  at  that  angle  from  linear  combinations  of  the 
corresponding  original  strip  functions,  where  these  linear  combinations  correspond  precisely  to  a 
1-dimensional  wavelet  transform  perpendicular  to  the  projection  direction.  This  transformation 
of  the  basis  ftmctions  is  shown  schematically  in  Figure  3  (which  corresponds  to  the  case  of  the 
rectangular,  Haar  wavelet).  The  original  strip  basis  functions  (rows  of  Tk)  are  illustrated  in  the  left 
half  of  the  figure,  while  the  corresponding  collection  of  multiscale  basis  functions  (rows  of  Tk)  are 
shown  in  the  right  half.  The  heavy  boimdaries  illustrate  the  support  extent  of  the  corresponding 
basis  element  while  the  or  ”  (together  with  shading)  notionaUy  indicate  the  sign  of  the 
function  over  this  region.  Note  that  the  number  of  basis  elements  in  the  original  (left  half)  and 
the  multiscale  (right  half)  framework  are  the  same,  as  they  must  be  since  the  multiscale  framework 
involves  an  orthonormal  chamge  of  basis.  We  may  naturally  group  the  multiscale  2-D  spatial  basis 
elements  into  a  hierarchy  of  scale  related  components  based  on  their  support  extent  or  spatial 
localization,  as  shown  in  the  figure.  The  basis  elements  defining  the  m-th  scale  in  such  a  group  axe 
obtained  from  the  rows  of  Tk  corresponding  to  (i.e.  scaled  by)  the  associated  wavelet  coefficients 
at  that  scale.  We  can  see  that  the  basis  functions  of  these  different  scale  components,  though  arising 
from  a  1-dimensional  multiscale  decomposition,  naturally  represent  behavior  of  the  2-dimensional 
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Figure  3:  Example  of  relationship  between  original  strip  basis  functions  contained  in  T*  (shown 
in  the  left  half  of  the  figure)  and  transformed  multiscale  basis  functions  of  7*.  (shown  in  the  right 
half  of  the  figure)  for  a  fixed  angle  k  corresponding  to  the  Haar  wavelet.  The  multiscale  basis 
functions  may  be  naturally  grouped  into  different  scale  components  based  on  their  spatial  extent 
(or,  equivalently,  their  relation  to  the  coefficients  in  ^jt),  as  shown. 

object  at  different  resolutions,  directly  corresponding  to  the  different  scale  components  contained 
in  the  transformed  vector  ik-  In  particular,  in  defining  the  overall  object  /,  the  multiscale  basis 
functions  at  scale  m  and  angle  k  are  weighted  by  the  corresponding  detail  component  The 

overall  object  is  then  represented  by  a  superposition  of  such  components  at  all  angles  k,  as  captured 
in  the  sum  in  (12). 

So  far  we  have  simply  transformed  the  representation  of  the  original  finest  scale  object  estimate 
/.  But  the  preceding  discussion  together  with  the  development  in  Section  2.3  suggests  how  to  use 
our  new  multiscale  decomposition  ^k  aJ^d  corresponding  basis  functions  7*.  to  obtain  a  multiscale 
decomposition  of  the  object  estimate  in  the  original  space.  Such  a  multiresolution  decomposition 
can  be  obtained  through  (12)  by  using  a  series  of  approximations  to  at  successively  finer  scales, 
thereby  inducing  a  series  of  corresponding  approximate  representations  of  the  object.  In  particulcir, 
we  define  the  m-th  scale  approximation  /I"*)  to  /  as; 

=  Y.T'HA{m)ik)  (13) 

*=1 

where  recall  that  the  m-th  scale  approximation  (A.(m)^fc)  is  obtained  by  zeroing  the  finer  scale 
components  in  the  vector  of  1-D  wavelet  transform  coefficients  of  as  discussed  in  Section  2.3. 
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Thus  the  approximation  uses  only  the  m  coarsest  scale  components  of  the  fuU  vector 
Similarly,  by  we  denote  the  additional  detail  required  to  go  from  the  object  approximation 

at  scale  m  to  that  at  scale  (m  +  1),  which  is  given  by: 

A 

(14) 

k=l 

where  recall  that  the  detail  vector  {D{m)  ^k)  is  obtained  by  zeroing  ail  but  the  corresponding  level 
of  detail  in  Combining  the  object  approximation  and  detail  definitions  (13)  and  (14) 
with  the  scale  recursive  relationship  (10)  we  see  that  the  object  itself  satisfies  the  following  scale 
recursive  relationship,  whereby  the  object  approximation  at  the  next  finer  scale  is  obtained  from 
the  approximation  at  the  current  (coarser)  scale  through  the  addition  of  the  incremental  detail  at 
this  scale,  just  as  for  the  1-D  case  treated  in  Section  2.3: 

f{m+l)  ^  ^ 

Note  that  our  multiscale  object  representation  given  in  (13)  and  corresponding  scale  recursive 
construction  (15)  is  induced  naturally  by  the  structure  of  the  individual  1-D  wavelet-based  mul¬ 
tiscale  decompositions  at  each  angle  k  and  is  not  simply  a  2-D  wavelet  tr2insform  of  the  original 
object  estimate  /.  In  other  words,  we  are  not  simply  relating  the  coefficients  of  a  2-D  multiscale 
decomposition  of  /  based  in  the  original  object  domain  to  those  of  a  1-D  decomposition  of  the  data 
in  the  projection  domain,  but  rather  we  are  allowing  a  multiscale  projection  domain  decomposition 
to  induce  a  corresponding,  and  thus  naturally  well  matched,  multiscale  object  representation.  In 
particular,  the  m-th  scale  approximation  of  the  object  is  created  as  a  linear  combination  of 
the  corresponding  m  coarsest  multiscale  basis  functions  (c.f.  Figure  3)  summed  over  all  singles  k 
(note  that  the  coefficients  finer  them  level  m  in  (A(m)  ^jt)  are  zero  and  use  the  object  definition 
(13)).  As  can  be  seen,  our  resulting  object  representation  lives  close  to  the  projection  domain  in 
which  data  is  gathered,  with  advantages  in  efficiency  as  we  will  see. 

Multiscale  Coefficient  Determination 

We  now  have  a  natural  multiscale  object  representation  framework  through  (13),  (14),  and  (15) 
that  is  similar  in  spirit  to  the  FBP  case  (3).  To  complete  the  process  and  create  multiscale  object 
estimates  from  data  we  must  find  the  multiscale  object  coefficients  (which  contain  all  the  infor¬ 
mation  we  need).  Further  we  desire  to  find  these  object  coefficients  directly  from  corresponding 
multiscale  tomographic  observations.  Aside  from  simply  being  am  evocative  notion  (e.g.  directly 
relating  scale-specific  data  features  to  corresponding  scale-specific  object  characteristics),  such  an 
approach  should  be  more  efficient,  in  that  we  would  expect  coarse  scale  object  characteristics  to 
be  most  strongly  affected  by  coarse  or  aggregate  data  behavior  and,  conversely,  fine  scale  object 
characteristics  to  depend  most  strongly  on  fine  scale  data  behavior.  Said  another  way,  we  would 
expect  the  relationship  between  such  multisczde  data  and  object  elements  to  be  nearly  diagonal, 
and  this  is  indeed  the  case. 

To  the  above  ends,  we  perform  a  wavelet-based  multiscale  change  of  basis  to  the  data  sequences 
yk,  similar  to  object  oriented  one  in  (11),  to  obtain  an  equivEdent  set  of  multiscale  observations: 

Vk  =  Wvk  (16) 

where,  recall,  W  is  a  matrix  taking  a  discrete  sequence  to  its  wavelet  tramsform.  We  may  now 
easily  obtaiin  our  desired  direct  relationship  between  the  multiscale  representation  of  the  data  at 
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angle  k  in  r/*  £ind  the  multiscsde  object  coefficients  at  the  same  angle  by  combining  the  two 
transformations  (11)  and  (16)  together  with  the  original  FBP  relation  (4)  to  obtain: 

ik  =  nr^k  (17) 

where  TZ.  =  RW  is  the  multiscale  data  filter,  corresponding  to  the  ramp  filter  R  of  the  usual  FBP 
case.  As  we  show  through  examples  later,  the  operator  R  is  compressed  by  the  wavelet  operator  so 
that  TZ  is  nearly  diagonal.  Further,  higher  compression  is  achieved  if  Daubediies  wavelets  jD„  with 
larger  n  are  used.  This  observation  is  consistent  with  the  observations  of  BeyUdn  et  al  [6],  since  R 
is  a  pseudo-differential  operator. 

The  Overall  Multiscale  Algorithm 

We  are  now  in  a  position  to  present  our  overall  multiscale  reconstruction  method.  By  comparing 
the  FBP  equations  (3)  and  (4)  to  the  corresponding  multiscale  equations  (12)  and  (17),  respectively, 
we  see  that  our  complete  multiscale  reconstruction  process  for  estimation  of  /  parallels  that  of  the 
standard  FBP  reconstruction,  in  that  identical  and  independent  processing  is  performed  on  the 
multiscale  data  sets  ijk  at  each  angle  to  obtain  the  corresponding  multiscsde  object  coefficients 
at  that  angle,  which  are  then  back-projected  along  corresponding  multiscale  basis  functions  7*.  and 
combined  to  obtain  the  final  object  estimate.  Thus  our  overall  procedure,  given  next,  is  no  more 
complex  than  the  standard  FBP  method. 

Algorithm  1  (Multiscale  Reconstruction) 

1.  For  a  given  choice  of  wavelet,  form  the  multiscale  filter  matrix  TZ  =  W  RW^  (the  multiscale 
counterpart  of  the  original  ramp  filter)  to  process  the  data  at  each  angle.  TZ  is  nearly  diagonal. 

2.  For  each  angle  k  perform  the  following: 

(a)  Find  the  multiscale  observations  tjk  by  taking  the  1-D  wavelet  transform  of  the  projection 
data  at  angle  k,  ijk  =  W yk- 

(b)  Calculate  the  multiscale  object  coefficient  set  ^k  =  'l^Vk  by  filtering  the  multiscale  obser¬ 
vations. 

(c)  Back-project  ^k  along  the  corresponding  multiscale  basis  functions  Tk,  T^^k- 

3.  Combine  the  object  contributions  from  the  individual  back-projections  at  each  angle  to  obtain 
the  overall  estimate,  J2k  ^k  ■ 

Beyond  simply  finding  a  finest  scale  object  estimate  as  described  in  Algorithm  1,  however, 
we  also  have  a  method  to  reconstruct  the  underlymg  object  at  multiple  resolutions  through  (13), 
(14)  amd  (15)  and  thus  for  easily  obtaining  information  about  the  object  at  midtiple  scales.  In 
particular,  if  an  approximation  at  scale  m  is  desired,  then  in  Algorithm  1  we  need  only  replace 

(,k  by  (A(m)  ^jt)  in  Step  2c  and  3.  In  particular,  this  simply  amounts  to  zeroing  detaul  components 
in  ^k  which  are  finer  than  scale  m.  Further,  if  instead  we  want  to  reconstruct  the  detail 
added  at  a  particular  scale,  we  need  only  replace  (,k  by  {D{rn)  ^k)  in  Step  2c  and  3  of  Algorithm  1. 
Similarly,  this  simply  amomits  to  zeroing  2ill  but  the  desired  scede  of  detail  in  £*.  Note  that 
such  intermediate  scale  information  about  /  can  even  be  efficiently  foimd  by  calculating  only  those 
elements  necessary  for  reconstructing  the  scale  of  interest  -  i.e.  aU  of  (,k  is  not  required.  For  example, 
if  8dl  that  is  required  is  a  coarse  estimate  of  the  object  and  not  the  fuU  reconstruction,  only  the 
coarsest  elements  of  are  required.  Conversely  if  only  fine  scale  features  are  to  be  reconstructed, 
then  only  the  finest  scale  detail  components  of  are  needed. 
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Figure  4:  Phantom  used  for  reconstruction  experiments.  The  phantom  is  256  X  256  and  projections 
axe  gathered  at  256  equally  spaced  angles  {Ns  =  256)  with  256  strips  per  angle  {Ng  =  256). 

Examples 

We  now  show  some  examples  of  om  multiscale  reconstruction  framework.  Figure  4  shows  the 
256  X  256  phantom  used  in  the  experiments  of  this  section.  Projection  data  were  collected  at  256 
equally  spaced  angles  {Ne  =  256)  with  256  strips  used  for  each  projection  {Ng  =  256).  First  we 
show  a  series  of  approximate  reconstructions  using  the  Daubechies  D3  wavelet  for  the  multiscale 
decomposition  W.  Figure  5  shows  the  various  scale  approximate  object  reconstructions  for 
the  entire  range  of  scales  m  —  1, . . . ,  8.  The  approximations  get  finer  from  left  to  right  and  top 
to  bottom  (so  that  the  upper  left  frcime  is  and  the  bottom  middle  frame  corresponds  to 
The  bottom  row,  right  shows  the  FBP  reconstruction  for  comparison.  Note  in  particular,  that 
the  finest  scale  approximation  is  identical  to  the  FBP  estimate  /.  The  intermediate  scale 
estimates  demonstrate  how  information  is  gathered  at  different  scales.  For  example,  in  the  scale 
3  reconstruction  (top  right  in  the  figure)  though  only  8  of  the  full  256  coefficient  elements  in 
the  vectors  are  being  used,  we  can  already  distinguish  separate  objects.  By  scale  4  (middle  row, 
left)  we  can  start  to  identify  the  septate  bright  regions  within  the  central  larger  object,  while  by 
scale  5  this  information  is  well  localized.  Even  at  this  comparatively  fine  scale  we  are  stiU  only 
using  about  12%  of  the  full  object  coefficient  set. 

In  Figure  6  we  show  the  corresponding  detail  components  for  the  same  phantom.  Again, 

the  additive  det2dl  becomes  finer  going  from  left  to  right  and  top  to  bottom  in  the  figure.  Notice 
that  the  fine  sceile,  edge  based,  features  of  the  phantom  are  clezirly  visible  in  the  and  A/*'®) 

reconstructions  (center  row,  middle  and  right  in  the  figure),  showing  that  structural  information  can 
be  obtained  from  these  detail  images  alone.  Recall  that  these  images  provide  the  added  information 
needed  in  going  from  the  object  approximation  at  one  scale  to  that  at  the  next  finer  scale  (as 
provided  in  Figure  5). 

As  we  discussed  earlier,  the  wavelet-based  multisczile  transformation  of  both  the  representation 
xit  and  data  also  serves  to  compress  the  ramp  filter  matrix  R  so  that  the  corresponding  mul¬ 
tiscale  filter  matrix  TZ  is  nejirly  diagonal.  As  we  argued  earlier,  this  reflects  the  fact  that  coarse 
scale  object  characteristics  are  most  strongly  affected  by  coarse  or  aggregate  data  behavior  and, 
conversely,  fine  scale  object  characteristics  tend  to  depend  most  strongly  on  fine  scale  data  behav¬ 
ior.  One  consequence  is  that  a  very  good  approximation  to  the  exact  reconstruction  procedure  of 
Algorithm  1  can  be  achieved  by  ignoring  the  off-diagonal  terms  of  TZ  in  (17).  These  off-diagonal 
terms  capture  both  intra  and  inter-scale  couplings.  Further,  this  approximation  to  the  exact  recon¬ 
struction  becomes  better  as  Daubechies  wavelets  D„  with  larger  n  are  used.  To  illustrate  this  point, 
in  Figure  7  we  show  complete  (finest  scale)  reconstructions  /  of  the  same  phantom  as  before,  based 
on  the  same  projection  data  but  using  a  diagonal  approximation  to  TZ  in  (17)  and  Algorithm  1  for 


Figiire  6:  The  detail  added  between  successive  scales  in  the  reconstructions  of  Figure  5.  First  row, 
left:  A/(0).  First  row,  middle:  A/(l).  First  row,  right:  A/(2).  Second  row,  left:  A/(3).  Second 
row,  middle:  A/(4).  Second  row,  right:  A/(5).  Third  row,  left:  A/(6).  Third  row,  middle:  A/(7). 


Figure  7:  Complete  finest  scale  multiscale  reconstructions  for  phantom  of  Figure  4  for  different 
approximate  filtering  operators.  The  left  three  frames  show  approximate  multiscale  reconstructions 
using  only  the  diagonal  elements  of  "R,  corresponding  to  different  choices  of  the  imderlying  wavelet: 
First  column  =  Haar.  Second  column  =  Dz-  Third  column  =  Ds-  For  comparison,  the  right-most 
frame  shows  an  equivalent  approximate  FBP  reconstruction  using  only  the  diagonal  elements  of  R, 
demonstrating  the  superiority  of  the  multiscale  based  approximations. 
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Figure  5:  Approximation  reconstructions  of  phantom  of  Figure  4  at  various  scales,  using  D3  wavelet. 
First  row,  left:  First  row,  middle:  First  row,  right:  Second  row,  left:  Second 

row,  middle:  Second  row,  right:  Third  row,  left:  Third  row,  middle:  The  third 

row,  right  shows  the  corresponding  FBP  reconstruction  /  for  comparison.  The  FBP  reconstruction 
is  the  same  as  since  this  is  the  complete  reconstruction. 

a  variety  of  choices  of  the  wavelet  defining  W.  For  the  reconstructions  we  use  only  the  diagonal 
elements  of  TZ  (which  account  for  0.0031%  of  aU  the  elements  for  this  case)  in  the  calculation  of 
effectively  setting  aU  off-diagonal  elements  to  zero.  Reconstructions  corresponding  to  Daubechies 
wavelets  £)„  with  increasing  n  (in  particular  Haar  or  Di ,  D3,  and  Dg)  are  shown  from  left  to  right 
in  the  figure.  It  can  be  seen  from  the  improvement  in  the  reconstructions  that  the  accuracy  of  the 
diagonal  approximation  becomes  better  as  wavelets  with  increasing  n  are  used  in  the  definition 
of  W.  In  particulzu:,  the  approximations  can  be  seen  to  compzire  very  favorably  with  the  standard 
FBP  reconstruction.  For  comparison  we  also  show  on  the  far  right  in  Figure  7  a  corresponding 
approximate  FBP  reconstruction  obtained  using  a  diagonal  approximation  to  the  original  ramp- 
filter  matrix  R  for  reconstruction.  It  can  be  seen  that  a  diagonal  approximation  in  the  multiscale 
domain  results  in  far  better  reconstructions  that  a  similar  approximation  in  the  original  domain, 
indicating  that  the  multiscale  transformation  of  data  and  coefficients  has  served  to  decouple  the 
resultant  quantities. 

In  sinnmary,  we  have  formulated  a  2-D  multiscale  object  reconstruction  method  in  terms  of 
approximation  and  detail  images.  This  method  is  derived  from  the  classical  FBP  method  and  thus 
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is  well  matched  to  reconstruction  from  projection  data.  The  associated  2-D  multiresolution  object 
representation  is  induced  by  a  1-D  wavelet-based  change  of  basis  to  the  original  FBP  projection 
space  object  coefficients.  While  the  resulting  representations  are  similar  in  spirit  to  a  direct  2-D 
multiresolution  decomposition  of  the  original  object,  in  that  approximations  are  produced  at  a  series 
of  scales  along  with  the  detail  necessary  to  proceed  from  one  such  approximation  to  the  next  finer 
one,  our  approach  does  not  correspond  to  such  a  direct  orthonormal  decomposition.  As  a  result  it 
is  fundamentally  different  from  previous  multiscale-related  work  in  tomography  (for  example,  [15]). 
In  these  approaches  such  a  direct  2-D  expansion  of  the  object  (i.e.  a  2-D  wavelet  transform)  is 
used  to  directly  define  the  approximation  and  detail  images,  the  coefficients  of  which  are  then 
calculated  from  the  projection  data.  In  contrast,  all  of  our  multiscale  quantities  inherently  “Uve” 
in  the  projection  domain.  As  a  result,  our  representation  is  closer  to  the  measurement  domain 
than  previous  multiscale  representations,  and  in  particular  implies  that  our  approach  is  no  more 
computation2illy  complex  than  FBP.  To  this  point  we  have  focused  on  noiseless  reconstructions. 
Next,  we  build  on  our  multiscale  reconstruction  method  to  obtain  a  fast  method  for  computing 
regularized  reconstructions  from  noisy  projections. 

4  Multiscale  Regularized  Reconstructions 

In  this  section  we  consider  the  estimation  of  an  object  /  from  noisy  projection  observations.  We  ex¬ 
tend  our  multiscale  reconstruction  method  presented  in  Section  3  to  obtain  statistically  regularized 
estimates  which  may  be  simply  and  efficiently  computed,  in  particular,  with  no  more  effort  them,  is 
required  for  the  standard  FBP  reconstruction.  This  regularized  solution  is  obtained  by  first  solving 
for  the  Maximum  Aposteriori  Probability  (MAP)  estimate  [29]  of  the  multiscale  object  coefficients, 
corresponding  to  a  certain  naturally  derived  multiscale  prior  model  and  then  back-projecting 
these  multiscale  coefficient  estimates  along  the  corresponding  multiscale  basis  ftmctions  as  before. 

The  presence  of  noise  in  projection  data  often  leads  to  reconstructions  by  standard  methods, 
such  as  FBP,  that  are  unacceptable  and  thus  require  some  form  of  regularization.  Traditionally,  two 
broad  approaches  have  been  used  in  the  generation  of  regularized  object  estimates  from  such  noisy 
projection  data.  Perhaps  the  simplest  approach  has  been  to  simply  roll  off  the  ramp  filter  used 
in  the  standzud  FBP  reconstruction  at  high  frequencies.  This  is  called  apodization  [7]  and  several 
different  windows  are  typically  used  for  this  pmrpose,  for  example  Hanning,  Hamming,  Parzen, 
Butterworth  etc.  [8].  The  assumption  is  that  most  the  object  energy  occurs  at  low  frequencies 
while  the  most  disturbing  noise-derived  artifacts  occur  at  high  frequency.  The  high  frequency  roll¬ 
off  thus  attenuates  these  components  at  the  expense  of  not  reconstructing  the  fine  scale  features 
in  the  object.  Since  the  overall  procedure  is  essentially  the  same  as  the  original  FBP  method, 
the  restdt  is  a  fast,  though  ad  hoc,  method  for  regularization.  The  other  traditional  approach  to 
regularizing  the  noisy  data  problem  is  statistically  based.  This  method  starts  with  a  statistical 
model  for  the  noisy  observations  based  on  (2): 

Vk  =  Tkf  4-  nfc  (18) 

where  n*.  is  tciken  as  an  additive  noise  vector  at  angle  k.  This  observation  model  is  then  coupled  with 
a  2-D  Markov  random  field  (MRF)  prior  model  [25,26]  for  /  to  yield  a  direct  MAP  estimate  of  the 
object  /.  While  statistically  based,  thus  allowing  the  systematic  inclusion  of  prior  information,  the 
2-D  spatiaUy-local  MRF  prior  models  used  for  the  object  generally  lead  to  optimization  problems 
that  are  extremely  computationedly  complex.  As  a  result,  these  methods  have  traditionally  not 
fotmd  favor  in  practical  apphcations. 
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In  contrast  to  the  above  two  techniques,  we  will  develop  a  multiscale  MAP  object  estimate  that, 
while  retaining  all  of  the  advantages  of  statistically-based  approaches,  is  obtained  with  the  same 
computational  complexity  as  the  FBP  reconstruction.  To  accomplish  this,  we  continue  to  work 
in  the  projection  domain,  as  the  FBP  method  does,  and  build  our  statistical  models  there,  rather 
than  in  the  original  object  domain.  As  in  Section  3,  we  then  allow  the  resulting  projection  domain 
coefficients  to  induce  a  2-D  object  representation  through  the  back-projection  and  summation 
operations.  To  this  end  we  start  with  an  observation  equation  relating  the  noisy  data  yk  to  the 
FBP  object  coefficients  Xk,  rather  than  the  corresponding  2-D  object  /  as  is  done  in  (18).  Such  a 
relationship  may  be  found  in  the  FBP  relationship  (4),  which  in  the  presence  of  noise  in  the  data 
becomes: 

yk  =  R~^Xk  +  nk,  ufe  ~  A/^(0,A„j)  (19) 

where,  recall  R  is  the  FBP  ramp  filter  operator^,  the  notation  z  ~  M'{m,A)  denotes  a  Gaussian 
distribution  of  meein  m  and  covariance  A  and  denotes  an  n  x  n  identity  matrix.  In  particular, 
we  assume  the  Anj  =  ^k^N,  >  he.  that  the  noise  is  uncorrelated  from  strip  to  strip  but  may  have 
different  noise  covariances  at  different  angles,  capturing  the  possibility  that  the  data  at  different 
projections  may  be  of  differing  quality  (e.g.  due  to  different  sensors  or  imaging  configurations). 
Further,  we  assume  that  the  noise  is  uncorrelated  from  eingle  to  angle,  so  that  n*  is  independent 
of  Uj,  k  ^  j.  This  model  of  independent  noise  in  the  projection  domain  is  well  justified  for  most 
tomographic  applications. 

As  in  Section  3,  for  purposes  of  estimation  we  desire  a  relationship  between  multiscale  repre¬ 
sentations  of  the  data,  object  coefficients,  and  noise.  Working  in  the  multiscade  transform  domain 
will  again  allow  us  to  obtain  induced  multiresolution  estimates  of  the  object.  Such  a  multiscale 
oriented  relationship  between  the  quantities  of  interest  can  be  found  by  combining  (19)  with  the 
multiresolution  changes  of  bases  (11)  and  (16)  based  on  W  (defined  in  Section  2.3)  to  obtain: 

Tjk  =  'R~^^k  +  Vk,  Vk  ~  A/'(0,  Ay^).  (20) 

where  Vk  =  Wuk  is  the  multiscale  trcinsformed  noise  vector  at  angle  k  with  A^.^  =  WAn^W^  = 
as  its  corresponding  covariance.  This  equation  relates  om  observed  noisy  multiscale  data  T)k 
to  our  desired  multiscale  object  coefficients  ^k  through  the  multiscale  filtering  operator  TZ.  Note  that 
the  assumption  of  uncorrelated  noise  from  angle  to  angle  and  strip  to  strip  in  the  original  projection 
domain  results  in  uncorrelated  noise  from  angle  to  angle  and  multiscale  strip  to  multiscale  strip  in 
the  multiscale  domain,  since  W  is  an  orthonormal  transformation. 

The  Multiscale  Prior  Model 

To  create  a  MAP  estimate  of  the  mrdtiscale  object  coefficients  $k,  we  will  combine  the  observation 
equation  (20)  with  a  prior  statistical  model  for  the  desired  unknown  multiscale  coefficient  vectors  $k- 
Multiresolution  object  estimates  and  the  detail  between  them  can  then  be  easily  obtained  by  using 
the  resulting  MAP  coefficient  estimates  ^k  ^.t  multiple  scales  in  the  multiscale  object  definitions 
(13)  and  (14),  as  was  done  previously  in  Section  3. 

We  base  our  prior  model  of  the  object  coefficients  directly  in  scale-space.  Such  scale-space 
based  prior  models  are  desirable  for  a  number  of  reasons,  e.g.  they  have  been  shown  to  lead  to 
extremely  efficient  scale-recursive  algorithms  [9, 13]  and  they  parsimoniously  capture  self-similar 

*Notc  that  (19)  assumes  that  R~^  exists.  For  the  case  where  R  represents  an  ideal  ramp  filter  this  wOl  indeed 
not  be  the  case,  as  this  operator  nulls  out  the  DC  component  of  a  signal.  For  filters  used  in  practice,  however,  this 
inverse  does  exist  and  the  expression  given  in  (19),  based  on  such  a  filter  is  well  defined.  Details  may  be  found  in 
Appendix  B. 
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behavior,  thus  providing  realistic  models  for  a  wide  range  of  natural  phenomenon.  In  particular, 
such  self-similar  models  can  be  obtained  by  choosing  the  detail  components  (i.e.  the  wavelet 
coefficients  at  each  scale)  as  independent,  A^(0,u^2”^)  reindom  variables  [14].  The  parameter 
p  determines  the  nature,  i.e.  the  texture,  of  the  resffiting  self-similar  process  while  cr^  controls 
the  oversill  magnitude.  This  model  says  that  the  variance  of  the  detail  added  in  going  from  the 
approximation  at  scale  m  to  the  approximation  at  scale  m  -f- 1  decreases  geometrically  with  scale. 
If  p  =  0  the  resulting  finest  level  representation  (the  elements  of  Xk)  correspond  to  samples  of  white 
noise  (i.e.  are  completely  uncorrelated),  while  as  p  increases  the  components  of  show  greater  long 
range  correlation.  SUch  self-similar  models  are  commonly  and  effectively  used  in  many  application 
areas  such  as  modeling  of  natural  terrain  and  other  textures,  biological  signals,  geophysiccJ  and 
economic  time  series,  etc.  [10-14]. 

In  addition  to  defining  the  scale  varying  probabilistic  structure  of  the  detail  components  of 
we  also  need  a  probabilistic  model  for  the  element  of  corresponding  to  the  coarsest  scale 
approximation  of  a;*,  i.e.  This  term  describes  the  DC  or  average  behavior  of  Xk,  of  which  we 
expect  to  have  little  prior  knowledge.  As  a  result  we  choose  this  element  as  A/’(0,Af),  where  the 
(scalar)  uncertainty  is  chosen  sufficiently  large  to  prevent  a  bias  in  our  estimate  of  the  average 
behavior  of  the  coefficients,  letting  it  be  determined  instead  by  the  data. 

In  summary,  we  use  a  prior  model  for  the  components  of  the  multiscale  coefficient  vectors 
which  is  defined  directly  in  scale-space  and  which  corresponds  to  a  self-similar,  fractal-like 
prior  model  for  the  corresponding  object  coefficients  Xk  •  lu  particular,  this  model  is  given  by 
~  ■A/^(0,  A()  with  $k  independent  from  angle  to  angle  and  where: 


A{  =  block  diag  Ja^'^  ,  A^^\  A]p,  Af  J 


(21) 


Agciin,  this  model  not  only  assumes  that  the  sets  of  multiscale  object  coefficients,  are  independent 
from  angle  to  angle  but  also  that  these  coefficients  are  independent  from  scale  to  scale,  that  they 
are  independent  and  identically  distributed  within  a  given  scale,  and  finally  that  their  variance 
decreases  geometrically  proceeding  from  coarse  to  fine  scales.  Obviously  other  choices  may  be  made 
for  the  statistics  for  the  multiscale  object  coefficients,  and  we  discuss  some  particularly  interesting 
possibilities  in  the  conclusions.  The  choice  we  have  made  in  (21)  while  simple,  is  well  adapted  to 
many  naturally  occurring  phenomenon.  In  addition,  since  the  observation  noise  power  is  uniform 
across  scales  or  frequencies,  the  geometrically  decreasing  variance  of  this  prior  model  implies  that 
the  projection  data  wUl  most  strongly  influence  the  reconstruction  of  coarse  scale  features  and 
the  prior  model  will  most  strongly  influence  the  reconstruction  of  fine  scale  features.  This  reflects 
our  belief  that  the  fine  scale  behavior  of  the  object  (corresponding  to  high  frequencies)  is  the 
most  likely  to  be  corrupted  by  noise.  Finally,  our  choice  of  prior  model  in  (21)  results  in  efficient 
processing  algorithms  for  the  solution  of  the  corresponding  MAP  estimate,  in  particular  with  no 
more  complexity  than  the  standard  FBP  reconstruction. 


The  Multiscale  MAP  Estimate 

We  are  now  in  a  position  to  present  our  overall  algorithm  for  computing  a  MAP  [29]  multisczile 
object  estimate  ^k-  Since  the  data  at  each  angle  rjk  and  the  corresponding  prior  model  for  ^k  are 
independent  from  angle  to  angle,  the  MAP  estimates  of  the  vectors  decouple.  In  particular,  the 
estimate  of  $k  at  each  angle,  based  on  the  observations  (20)  and  the  prior  model  (21)  is  given  by: 

ik  =  [Aj'  -H  7^-^A-'  Vk  =  nvk  (22) 
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where  the  regularized  multiscale  filter  operator  TZ  is  defined  in  the  obvious  way.  This  regularized 
filtering  matrix  is  exactly  analogous  to  the  unregularized  filtering  operator  7^  of  (17)  for  the  noise 
free  case.  In  this  regularized  case,  however,  TZ  now  also  depends  on  both  the  noise  model 
and  the  prior  object  model  A^.  If  the  noise  variance  is  low  relative  to  the  uncertainty  in  the  prior 
model  (so  is  large)  then  TZ  will  approach  TZ  and  the  estimate  will  tend  toward  the  standard 
unregularized  one.  Conversely,  as  the  noise  increases,  TZ  will  depend  to  a  greater  extent  on  the 
prior  model  term  Aj  and  the  solution  will  be  more  regularized  or  smoothed. 

Finally,  as  in  the  noise-less  case,  the  resulting  object  estimate  f  is  then  obtained  by  back- 
projecting  the  estimated  multiscale  object  coefficients  along  the  corresponding  multisc2de  basis 
functions  Tk  and  combining  the  result.  The  overall  structure  of  this  regularized  reconstruction 
parallels  that  of  the  original  FBP  method,  and  therefore  is  of  the  same  computational  complexity 
as  FBP.  In  summeiry,  om  overall,  efficient  regularized  multiscale  estimation  algorithm  is  given  by 
the  following  procedure,  which  paredlels  our  unregularized  multiscale  reconstruction  algorithm: 

Algorithm  2  (Regularized  Multiscale  Reconstruction) 

1.  Find  the  regularized  multiscale  filter  matrix  TZ  (the  multiscale  regularized  counterpart  of  the 

original  ramp  filter)  by  doing  the  following: 

(a)  For  a  given  choice  of  wavelet,  form  the  unregularized  multiscale  filter  matrix  TZ  = 
WKW"^  as  before. 

(b)  Choose  the  model  parameters  Xk  specifying  the  variances  of  the  observation  noise  pro¬ 
cesses  and  thus  defining  c.f.  (20). 

(c)  Choose  the  multiscale  prior  model  parameters  p  and  A(  specifying  the  magnitude  and 

texture  of  the  model  and  the  uncertainty  in  its  average  value,  respectively,  and  generate 
the  prior  covariance  matrix  A(  through  (21). 

(d)  FormTZ=  [a^^^  +TZ-^A-^TZ-^Y^ 

2.  For  each  angle  k  perform  the  following: 

(a)  Find  the  multiscale  observations  Tjk  by  taking  the  1-D  wavelet  transform  of  the  projection 
data  at  angle  k,  7}k  =  W y^. 

(b)  Calculate  the  regularized  multiscale  object  coefficient  set  =TZ‘qk  by  filtering  the  mul¬ 
tiscale  observations. 

(c)  Back-project  along  the  corresponding  multiscale  basis  functions  Tk,  T^^k- 

3.  Combine  the  regularized  object  contributions  from  the  individual  back-projections  at  each  angle 

to  obtain  the  overall  regularized  object  estimate,  Ylk'^ 

As  before,  we  may  also  easily  obtain  regularized  reconstructions  of  the  object  at  multiple  resolutions 
by  using  (13)  and  (14)  together  with  the  MAP  coefficient  estimates  £,k-  In  p2irticular,  to  obtmn 
the  approximation  at  scale  m  then  we  need  only  replace  ik  by  (corresponding  to 

simply  zeroing  some  of  the  terms  in  ^k)  in  Step  2c  and  3.  Similarly,  the  corresponding  object  detail 
components  at  scale  m  may  be  obtained  by  using  {D{m)^k)  in  place  of  ^k  in  these  steps. 

While  Algorithm  2  is  already  extremely  efficient,  m  that  2-D  multiscale  regularized  object  es¬ 
timates  are  generated  with  no  more  complexity  than  is  needed  for  the  standard  FBP  method, 
additional  gains  may  be  obtained  by  exploiting  the  ability  of  the  wavelet  transform  operator  W 
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to  compress  the  FBP  filtering  operator  R.  Recall,  in  particular,  that  the  (unregularized)  multi¬ 
scale  filtering  matrix  TZ  =  W RW^  is  nearly  diagonal,  with  this  approximation  becoming  better  as 
Daubechies  wavelets  with  increasing  n  are  used  in  the  specification  of  W.  Based  on  our  assump¬ 
tions,  the  matrices  and  specifying  the  prior  model  and  observation  covariances  respectively, 
are  already  diagonal.  If  in  addition  TZr^  were  also  a  diagonal  matrix,  then  from  (22)  we  see  that  TZ 
itself  would  be  diagonal,  with  the  result  that  the  “filtering”  in  Step  2b  of  Algorithm  2  would  simply 
become  point  by  point  scaling  of  the  data.  To  this  end  we  will  assume  that  the  wavelet  transform 
W  truly  diagonalizes  R  by  effectively  ignoring  the  small,  off-diagonal  elements  in  TZ~^.  That  is,  we 
assume^; 

7^"^  w  diag(ri,r2,...,rjvJ  (23) 

where  Vi  are  the  diagonal  elements  of  1Z~^ .  Now  let  us  represent  the  diagonal  prior  model  covariance 
matrix  as  A^  =  diag[pi,p2)  •  •  •  jPiV,])  and  recall  that  =  Xk^N,-  Using  these  quantities  together 
with  our  approximation  to  in  the  specification  of  the  estimate  (22)  yields  an  approximate 
expression  for  ^k- 


$k  w  diag 


ri 


r2 


I*iV, 


/i  +  i^klpiY  rl  -f-  (A*/p2)’  ■■■’  +  {XkIPN.), 


Vk  =  '^Vk 


(24) 


where  the  approximate  MAP  filtering  matrix  1Z  is  defined  in  the  obvious  way.  Our  experience  is 
that  when  W  is  defined  using  Daubechies  wavelets  of  order  3  or  higher  (i.e.  using  D3,  D4,...),  the 
estimates  obtained  using  'R  in  place  of  the  exact  regularized  filter  TZ  in  Algorithm  2  are  visually 
indistinguishable  from  the  exact  estimates  where  TZ~^  is  not  assumed  to  be  diagonal.  Indeed,  it  is 
actuaUy  this  approximate  filtering  operator  TZ  that  we  use  to  generate  the  example  reconstructions 
we  show  next. 

Before  proceeding,  however,  let  us  examine  our  MAP  reguleirized  filtering  operator  71  in  more 
detail  to  understand  how  our  multiscale  MAP  estimation  procedure  relates  both  to  the  standcird 
FBP  method  and  the  ad  hoc  regularization  obtained  through  apodization  of  the  FBP  filter.  The 
MAP  estimates  ^k  induce  corresponding  estimates  Xk  of  the  original  object  coefficients  Xk  through 
the  change  of  basis  (11)  and,  similarly,  rjk  and  y*  are  related  through  (16).  Thus,  the  multiscale 
MAP  estimation  operation  specified  by  (22)  imposes  a  corresponding  relationship  between  the 
original  finest  scale  quantities  Xk  and  y*,  given  by: 

Xk={w^7ZW)yk  =  R^yk  (25) 

where  the  effective  multiscale  MAP  regularized  filtering  matrix  Reff  is  defined  in  the  obvious  way. 
The  effect  of  this  MAP  regularized  filter  can  now  be  compared  to  the  standard  FBP  or  apodized 
ones.  The  behavior  of  the  matrix  operator  JKeff  catn  be  most  easily  imderstood  by  examining  its 
corresponding  frequency  domain  characteristics.  To  this  end,  in  Figure  8  we  plot  the  magnitude 
of  the  Fourier  transform  of  the  central  row  of  effective  regularized  matrix  Reg  corresponding  to  a 
variety  of  choices  of  the  model  or  regularization  parameters  A*,  (the  noise  variance)  and  p  (the  decay 
rate  across  scales  of  the  added  detail  variance)  for  fixed  =  1  (overall  prior  model  amplitude) 
and  A(  =  1  (prior  model  DC  variance).  We  also  plot,  with  heavy  lines,  the  magnitude  of  the 
Fourier  transform  of  the  corresponding  central  row  of  the  standard  FBP  ramp  filter  matrix  R  for 
comparison.  Prom  Figure  8,  we  can  see  that  in  the  multiscale  MAP  framework  regularization  is 

®One  can  imagine  another  level  of  approximation  where  we  set  the  off-diagonal  elements  of  'll  itself  to  zero  prior 
to  inversion  rather  than  those  of  7l~^.  This  further  approximation  results  in  reconstructions  which  are  visually  very 
similar  to  what  we  obtain  here. 
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Figtire  8:  The  Fourier  transform  of  the  central  row  of  iZeff  for  different  values  of  regularization 
parameters  p  and  A*,,  illustrating  the  effect  of  the  multiscede  regulctrizing  filter  in  the  frequency 
domain.  In  each  of  the  subplots,  the  V-shaped  heavy  line  corresponds  to  the  stcmdard  FBP  ramp 
filter  and  the  four  curves  from  top  to  bottom  correspond  to  p  =  0.5  (solid  line),  1.0  (dashed  Kne), 
1.5  (dashdot  line)  smd  2.0  (dotted  line)  respectively  (in  some  subplots  some  of  the  lines  overlap). 
In  all  cases  we  fixed  =  1  (the  overall  prior  model  amplitude)  and  =  1  (the  prior  model  DC 
variance). 

basically  achieved  by  rolling  off  the  ramp  filter  at  high  frequencies,  the  same  principle  as  used  in  the 
ad  hoc,  apodization  regularized  FBP  reconstructions.  We  also  see  that  decreasing  the  observation 
noise  variance  Afe  for  a  fixed  prior  model  structure  p,  or  conversely,  increasing  the  variance  of  the 
detail  added  in  proceeding  from  coarse  to  fine  scales  in  the  prior  model  (i.e.  decreasing  p)  for  a 
fixed  observation  noise  variance  Afc,  leads  to  decreased  regularization  as  reflected  in  decreased  high 
frequency  attenuation.  This  behavior  is  reasonable,  in  that  in  the  first  case,  the  data  becomes  less 
noisy  while  in  the  second  the  uncertainty  in  the  prior  model  becomes  larger.  In  both  these  cases 
one  would  Weint  to  put  more  reliance  on  the  data  (i.e.  less  regularization). 

In  summary  then,  our  multiscale  based  regiilarization  approach,  though  derived  from  statistical 
considerations  and  possessing  aU  the  advantages  of  such  methods  (e.g.  the  ability  to  incorporate 
prior  knowledge  in  a  rationrd  way,  the  ability  to  do  performance  amalysis  and  understand  the 
relative  importance  of  various  sources  of  imcertainty,  etc.),  obtains  residts  at  no  greater  (and  in 
some  cases  with  substantially  less)  computational  complexity  than  standard  unregularized  or  ad 
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hoc  approaches.  In  addition,  we  obtain,  essentially  for  free,  estimates  at  multiple  resolutions  and 
thus  the  ability  to  extract  information  from  data  at  multiple  scales. 

Examples 

Next  we  show  some  examples  of  reconstructions  using  our  multiscale  methods  in  the  presence  of 
noise.  The  scime  256  X  256  phantom  shown  in  Figure  4  was  used  for  all  experiments.  In  each 
case  projection  data  for  the  phantom  were  ageiin  generated  at  Ns  =  256  equally  spaced  angles 
with  Ng  =  256  strips  in  each  projection.  These  noise-free  values  were  then  corrupted  through  the 
addition  of  independent,  zero-mean  Gaussian  noise  to  yield  our  observations.  The  variance  A„  of 
this  additive  noise  depended  on  the  experiment  aind  was  chosen  to  yield  an  equivalent  signal-to-noise 
ratio  (SNR)  of  the  resulting  observations,  defined  as: 

SNR(dB)  =  101og^^i|^|®  (26) 

XnNsNg 

where,  recall,  Tkf  is  the  noise-free  projection  data  at  angle  k.  Finally,  in  all  multiscale  reconstruc¬ 
tions  we  show  here  the  Daubechies  £>3  wavelet  was  used  in  the  definition  of  W  for  the  reconstruction. 

The  first  example,  shown  in  Figure  9,  demonstrates  reconstruction  from  noisy  data  using  the 
unregularized  multiscale  approach  of  Section  3.  The  varizmce  An  of  the  added  noise  was  chosen  to 
yield  a  SNR  of  5  dB.  This  figure  shows  the  various  scale  approximate  object  reconstructions 
corresponding  to  the  unregularized  Algorithm  1  for  the  complete  range  of  scales  m  =  1, . . .,  8.  As 
before,  the  approximations  become  finer  from  left  to  right  and  top  to  bottom  (so  that  the  upper  left 
frame  is  and  the  bottom  middle  frame  corresponds  to  /'^*^).  The  bottom  right  frame  shows  the 
standard  FBP  reconstruction  based  on  the  noisy  data.  Since  corresponds  to  the  unregularized 
complete  finest  scale  reconstruction  it  is  also  the  same  as  the  standard  FBP  reconstruction  based 
on  the  noisy  data  for  this  case.  The  figure  illustrates  the  resolution-accuracy  tradeoff  inherently 
captured  in  the  multiscale  framework  and  confirms  the  point  that  even  in  the  unregularized  case, 
information  from  noisy  observations  can  be  focused  by  stopping  the  reconstruction  at  a  coarse  scale, 
for  example  scale  5  (center  row,  middle  in  the  figure).  The  finer  scale  detail  contributions 
m  >  5  are  evidently  mainly  noise  which  obscure  the  object  features.  In  particular,  in  the  finest 
scale  reconstruction  (i.e.  the  standcird  FBP  reconstruction)  the  object  is  zdmost  completely  lost  in 
the  noise. 

Next  we  show  estimates  generated  by  our  multiscale  MAP  regularized  estimation  method  dis¬ 
cussed  in  this  section.  Figure  10  shows  the  various  scale  approximate  object  reconstructions  /f"*) 
corresponding  to  our  multiscale  MAP  estimate  of  using  noisy  data  with  same  SNR  (i.e.  SNR  = 
5  dB)  as  in  Figure  9.  The  MAP  estimate  was  generated  using  the  extremely  efficient  approximate 
expression  (24),  which,  for  the  Daubechies  Dz  wavelet  we  are  using,  was  indistinguishable  from  the 
corresponding  estimate  based  on  the  exact  expression  (22).  Again  the  approximations  become  finer 
from  left  to  right  and  top  to  bottom  in  the  figure.  For  these  reconstructions  we  chose  the  modeled 
observation  noise  vmance  as  A*  =  5.5  X  10®.  For  the  statistical  model  parameters  of  the  prior,  the 
decay  rate  across  scale  of  the  added  detail  vEiriance  was  chosen  as  p  =  1.5,  the  overall  magnitude 
of  the  prior  was  set  to  cr^  =  11,  and  the  variance  of  the  prior  model  average  value  was  Aj  =  1. 
The  efiect  of  the  regidarization  can  be  readily  seen  in  its  ability  to  suppress  noise  in  the  finest 
scale  reconstruction.  For  comparison,  the  standard  FBP  reconstruction  for  this  case  is  given  on  the 
bottom  row,  right  in  Figure  10.  In  addition,  the  multiscale  nature  of  the  information  focusing  can 
be  seen  in  the  scale  evolution  of  the  estimates.  In  particulrir,  there  appears  to  be  little  difference 
between  scale  5  and  finer  scale  estimates  in  the  figure,  suggesting  that  little  additional  information 
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Figme  9:  Reconstructions  of  phantom  of  Figure  4  from  5  dB  SNR  projection  data  based  on  unreg¬ 
ularized  Algorithm  1  using  Dz  wavelet.  Reconstructions  are  shown  at  various  scades  demonstrating 
the  smoothing  eifect  that  can  be  achieved.  First  row,  left:  First  row,  middle:  First 

row,  right:  Second  row,  left:  Second  row,  middle:  Second  row,  right:  Third 

row,  left:  Third  row,  middle:  The  standaird  FBP  is  shown  in  the  third  row,  right  for 

comparison.  The  FBP  reconstruction  is  the  same  as  since  this  is  the  complete  unregularized 
reconstruction. 

is  being  obtained  in  proceeding  to  such  finer  scales,  that  the  additional  degrees  of  freedom  being 
added  at  such  finer  scales  are  not  really  being  supported  by  the  data,  and  thus  that  we  should 
stop  the  reconstruction  at  this  coarser  scade.  Further,  estimates  at  scale  5  and  coairser  appeair  quite 
similar  to  the  corresponding  unregnlarized  estimates  in  Figure  9,  showing  that  these  coarser  scaile 
estimates  are  dominated  by  the  data  and  are  not  very  dependent  on  the  prior  model  at  this  point 
anyway. 

Finadly,  in  Figure  11,  we  show  a  series  of  finest  scade  multiscale  MAP  regtdarized  reconstructions, 
corresponding  to  different  choices  of  the  prior  model  texture  as  determined  by  the  parameter  p. 
The  same  phaintom  as  before  is  used,  but  we  use  observations  with  a  SNR  of  —10  dB  (much 
worse  thain  used  above).  The  standaurd  FBP  reconstruction  is  shown  for  compairison  in  the  far 
right  image  of  the  figure.  The  object  is  completely  lost  in  the  FBP  reconstruction  at  this  extreme 
level  of  noise.  The  MAP  reconstructions  aue  shown  in  the  first  three  frames  of  the  figure,  with  a 
smoother,  more  correlated  prior  model  being  used  as  we  proceed  from  left  to  right.  The  specific 
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Figttre  10:  Multiscale  MAP  regulzirized  reconstructions  at  various  scales  of  phemtom  of  Figure  4 
from  5  dB  SNR  projection  data  using  Dz  wavelet.  The  values  of  the  statistical  model  parameters 
used  are  Xk  =  5.5  X  10®,  p  =  1.5,  cr^  =  11,  Af  =  1.  First  row,  left:  First  row,  middle: 

First  row,  right:  Second  row,  left:  Second  row,  middle:  Second  row,  right: 

Third  row,  left:  Third  row,  middle:  For  comparison,  the  stemdard  FBP  reconstruction 

for  this  case  is  given  iu  the  third  row,  right.  The  improved  ability  of  the  regularized  reconstructions 
to  extract  information  is  demonstrated. 

multiscale  MAP  model  parameters  were  chosen  as  follows.  The  observation  noise  variemce  was 
chosen  as  A*  =  1.7  X  10^.  The  overall  prior  model  magnitude  was  set  to  cr*  =  17  while  the  prior 
model  DC  vEiriance  was  set  to  Af  =  1.  The  prior  model  texture  parameter  p  took  on  the  values 
{0.5, 1.0, 1.5}.  The  increased  smoothness  in  the  prior  can  be  seen  to  be  reflected  in  increased 
smoothness  of  the  corresponding  estimates.  Note  also  the  ability  of  the  algorithm  to  pidl  out  at 
least  the  global  object  features  in  the  presence  of  this  substanticd  amount  of  noise.  Again,  the  more 
highly  smoothed  reconstructions  (corresponding  to  higher  values  of  p)  appear  quite  similar  to  the 
coarser  level,  unregularized  reconstructions  shown  previously,  showing  that  we  are  really  accessing 
the  cosuse  level  information  in  the  data. 


25 


Figure  11:  Multiscale  MAP  regularized  reconstructions  of  the  phantom  of  Figure  4  at  the  finest  scale 
from  —10  dB  SNR  observations  for  different  choices  of  prior  model  texture,  p,  with  Afe  =  1.7  X  10^, 

=  17,  and  Aj  =  1,  are  shovm  in  the  first  three  frames:  Frame  p  —  0.5.  Frame  2:  p  =  1.0. 
Frame  3:  p  =  1.5.  For  comparison  the  standard  FBP  reconstruction  is  shovra  in  the  last  frame  on 
the  far  right. 

5  Conclusions 

In  this  paper  we  have  developed  a  wavelet-based  multiscale  tomographic  reconstruction  technique 
which  is  different  from  other  multiscale  techniques  in  the  following  respects.  First,  our  2-D  mul¬ 
tiscale  object  representation  is  naturally  induced  by  expanding  the  FBP  coefficients,  and  hence 
basis  functions  (i.e.  strips),  in  a  1-D  wavelet  basis.  This  is  in  contrast  to  other  multiscale  recon¬ 
struction  techniques  which  begin  with  a  2-D  object  representation  obtained  from  a  full  2-D  wavelet 
decomposition  of  the  object  space.  These  techniques  must  subsequently  relate  the  inherently  1-D 
projection  data  to  these  fundamentally  2-D  object  coefficients.  In  contrast,  the  multiscale  repre¬ 
sentation  resulting  from  our  approach,  arising  as  it  does  from  the  projection  strips  themselves,  is 
much  closer  to  the  measurement  domain.  The  result  is  a  highly  efficient  method  to  compute  our 
multiscale  object  coeflScients,  m  particular,  no  more  complex  than  the  widely  used  standard  FBP 
operation.  Yet,  rmlike  the  FBP  method,  our  multiscale  reconstructions  also  provide  a  framework 
for  the  extraction  and  presentation  of  information  at  multiple  resolutions  from  data.  Further,  our 
resulting  multiscale  relationships  between  data  and  object  allow  extremely  simple  approximations 
to  be  made  to  our  exact  relationships  with  virtually  no  loss  in  resulting  image  quality,  thus  further 
improving  the  potential  efficiency  of  our  approach.  Such  approximations  are  not  possible  with  the 
standard  FBP  method,  as  they  result  in  severe  artifacts. 

In  addition,  based  on  this  wavelet-based  multiscale  frcimework,  we  have  proposed  a  statistically- 
based  multiresolution  MAP  estimation  algorithm.  This  method  provides  statistically  regularized 
reconstructions  from  noisy  data,  and  does  so  at  mtdtiple  resolutions,  at  no  more  effort  than  is 
required  for  the  standard  FBP  method.  This  approach,  based  on  the  construction  of  prior  mod¬ 
els  directly  in  scale-space,  allows  for  the  inclusion  of  natural,  self-similar  prior  models  into  the 
reconstruction  process.  In  contrast,  conventional  statisticadly-based  regularization  methods,  utiliz¬ 
ing  MRF-type  prior  models  constructed  directly  in  (finest  scale)  object  space,  lead  to  extremely 
complex  and  taxing  optimization  problems.  The  result  has  typically  been  that  such  statistically 
motivated  methods  have  been  shTinned  in  practice  in  favor  of  fast,  though  ad  hoc,  approaches.  Our 
results  provide  a  bridge  between  these  two  extremes.  Further,  in  providing  estimates  at  multiple 
resolutions,  our  results  provide  tools  for  the  assessment  of  the  resolution  versus  accuracy  tradeoff, 
wherein  we  expect  coarser  scale  features  of  data  to  be  more  accurately  determined  than  finer  scale 
ones.  Though  we  did  not  exploit  this  ability  in  the  present  paper,  our  formulation  also  allowed 
for  the  possibility  of  combining  data  from  projections  of  fundamentally  different  quality,  through 


the  specification  of  different  noise  variances  A*,  at  different  angles.  The  resulting  estimates  do  not 
correspond  to  a  simple  FBP  or  even  roUed  off  FBP  reconstruction,  yet  are  easily  obtain  in  our  frame¬ 
work.  Finally,  as  before,  our  multiscale  MAP  approach  leads  to  algorithms  which  are  amenable  to 
an  additional  level  of  approximation,  with  resulting  improved  efficiency,  again  at  virtually  no  loss 
in  corresponding  reconstruction  quality. 

Even  though  this  paper  concentrates  on  the  complete-data  case,  our  multiscale  reconstruction 
method  has  the  potential  of  overcoming  this  limitation  and  dealing  with  limited  or  missing  data 
problems.  First  note  that  our  ability  to  combine  projections  of  different  quality  already  does  this 
to  some  extent,  in  that  some  projections  can  be  down  weighted.  More  important,  however,  is  the 
structure  of  the  prior  model  covariance.  Specifically,  while  the  prior  model  (21)  assumes  that  the 
multiscale  object  coefficients,  are  independent  from  auigle  to  angle,  we  would  intuitively  expect 
the  coarse  scale  object  coefficients  at  different  projection  angles  to  actually  be  highly  correlated 
with  each  other,  and  further  for  this  correlation  to  decrease  at  finer  scales.  Such  a  correlation 
structure  across  projection  angles  would  help  us  estimate  at  least  the  coarse  scale  object  coefficients 
to  a  good  accuracy  even  if  the  projection  data  at  certain  angles  are  missing.  The  current  angular 
independence  of  the  coefficients  corresponds  to  ein  overall  covariance  structure  for  these  variables 
(i.e.  the  vector  of  all  £fc’s)  which  is  block  diagonal,  and  it  is  this  block  diagonality  that  is  partially 
responsible  for  the  extreme  efficiency  of  our  current  technique.  Using  the  proposed,  more  correlated 
prior  covariance  structure  would  correspond  in  this  paradigm  to  the  addition  of  off-diagonal  terms. 
At  first,  such  a  proposal  would  seem  to  make  things  dramatically  worse  from  a  computational 
perspective,  since  we  must  now  contend  with  what  corresponds  to  the  inversion  of  a  full  rather 
than  a  block-diagonal  matrix.  All  is  not  lost,  however,  for  at  least  two  reasons,  both  related  to 
the  fact  that  we  are  building  our  prior  models  directly  in  scale  space.  First,  the  addition  of  only 
coarse  scale  correlations  may  be  sufficient  to  regularize  a  given  problem,  with  the  result  that  only 
a  few,  low  dimensional,  off  diagonal  elements  need  be  added  to  the  prior  model  covariance  (recall, 
at  coarser  scales  there  are  far  fewer  model  elements  -  e.g.  at  the  coarsest  sczde  there  is  only  one 
per  angle).  These  few  additional  coarse  scale  terms  could  then  be  aggregated  into  a  slightly  larger 
corresponding  covariance  block,  returning  us  to  the  block  diagonal  case,  but  with  one  block  slightly 
larger  than  the  rest.  More  significantly  perhaps,  however,  is  that  recent  research  has  demonstrated 
that  certain  scale-based  prior  models  (which  correspond  to  tree  structures),  though  corresponding  to 
highly  correlated  fields,  can  lead  to  extremely  efficient  scale-rectirsive  estimation  algorithms  [9,13]. 
We  excimine  such  variations  of  prior  covariance  structure,  including  the  combination  of  projections 
of  differing  quadity  in  a  future  paper. 


27 


A  Summary  of  Notation 


Table  1  summarizes  the  notation  used  in  this  paper. 


Notation 

Explanation 

Ns 

Number  of  angular  projections. 

N. 

Number  of  strip  integrals  in  each  cingular  projection. 

k,l 

FBP  quantities  are  indexed  by  k,i. 
k  is  the  angle  of  projection,  A:  =  1, . . . ,  Ns. 

1  is  the  strip  within  the  angular  projection,  1  =  1, . .  .,Ng. 

f 

The  discretized  object  defined  on  a  iV,  x  N,  square  grid. 

f 

The  reconstructed  object. 

Vk 

Projection  data  set  at  angle  k,  k  =  yk  =  [yfc(l) ^*(2)  . . .  yk{N,)]'^. 

FBP  object  coefficient  set  at  Eingle  k,  k  =  1,. . Ns.  Xk  =  [**(!)  ®fe(2)  . . . 

Tk 

The  discrete  strip  projection  operator  at  angle  k,  yk  =  Tkf- 

n 

The  back-projection  operator  at  angle  k,  f  =  Ylk'^k^k- 

R 

The  matrix  representing  the  FBP  ramp  filtering  operation,  Xk  =  Ryk- 

W 

The  matrix  realization  of  the  discrete  1-D  wavelet  transform  operation. 

Dn 

Daubechies  wavelet,  where  n  indicates  the  support  size  of  the  corresponding  filter. 

m 

The  multiscale  quantities  are  indexed  in  scale  by  m; 

Scales  become  finer  with  increasing  m. 

Vk 

The  1-D  wavelet  transform  of  yk,  rjk  =  Wyk- 

6 

The  1-D  wavelet  transform  of  Xk,  ik  =  ^  ^ 

tC"*) 

Detail  vector  of  Xk  at  scale  m  containing  the  wavelet  transform  coefficients  of  Xk 
at  scale  m. 

Coarsest  scale  approximation  of  Xk- 

Tk 

The  multiscale  projection  operator  at  angle  k,  7]k  =  Tkf,  Tk  =  WTk. 

The  multiscale  back-projection  operator  at  angle  k,  f  =  Ylk  '^^k- 

f{^) 

The  approximate  object  reconstruction  at  scale  m. 

The  object  detail  reconstruction  at  scale  m,  =  /f'")  -f 

n 

The  multiscale  filter,  ^k  =  N-rik,  71  =  WRW^ . 

p,  h 

The  MAP  prior  model  parameters; 
p  is  the  decay  rate  of  prior  variance  across  scales. 

(T^  is  the  overall  prior  model  amplitude, 
is  the  prior  model  DC  variance. 

nk 

The  noise  vector  at  angle  k,  yk  =  R~^Xk  -f  nfel 

Uk  ~  A/’(0,  A„j),  A„j  =  XkiN.- 

Vk 

The  1-D  wavelet  transform  of  Uk,  Vk  =  Wuk,  Vk  ~  A/’(0,  A„j),  A^.^  =  Xk^N,', 

Vk  =  7l~^$k  +  Vk. 

n 

The  regularized  multiscale  filter. 

n 

An  approximation  to  above  obtained  by  ignoring  off-diagonal  elements  in  71. 

Table  1:  Notation  used  in  this  paper. 
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B  Details  about  the  formation  of  FBP  ramp-filter  matrix  R 

In  this  work  we  take  the  matrix  R  to  represent  the  practical  FBP  filtering  operator.  In  the  ideal 
case,  this  FBP  ramp  filter  operation  is  given  by: 

HnFn  (27) 

where  Fjv  is  a.  N  x  N  matrix  representing  the  1-D  Fourier  transform  operation  on  a  sequence  of 
length  N,  and  Hjf  is  a,  N  x  N  diagonal  matrix  containing  ideal  high-pass  “ramp”  filter  coefficients 
for  a  length  N  sequence.  The  matrix  Hjf  has  a  diagonal  entry  of  0  since  it  gives  zero  weight  to  the 
frequency  cell  centered  aroimd  0.  Thus  the  ideal  raunp  filter  coeflSicient  matrix  Hn,  and  hence  the 
matrix  (27),  is  not  invertible.  In  practice,  however,  to  avoid  dishing  (i.e.  interperiod  interference) 
and  DC  artifacts,  a  filtering  operator  R  is  used  which  is  constructed  according  to  [l]: 

R=SF^^H2nF2nS^  (28) 

where  F^n  is  a  2N  X  2N  matrix  representing  the  1-D  Fourier  transform  operation  on  a  sequence 
of  length  2N,  Huf  is  a  2N  X  2N  diagonal  matrix  containing  the  corresponding  ideal  ramp  filter 
coefficients,  and  the  N  X  2N  zero-padding  matrix  S  is  given  by 

5  =  [  0  In  0  ]  .  (29) 

If  we  define  H  to  be  the  equivalent  N  X  N  practical  ramp  filter  coefficient  matrix  such  that: 

R  =  F^^HFn  =  S  F^^  H2n  F2N  (30) 

then  H  can  be  seen  to  be  given  by: 

H  =  FnS  F-^  H2n  F2N  F^^  .  (31) 

One  can  show  that  H  is  a  diagonal  matrix  and  the  diagonal  elements  of  H  are  almost  identiccd  to 
Hff  except  that  H  gives  small  positive  weighting  to  the  frequency  cells  centered  around  0  (refer 
to  [1]  for  a  plot  of  the  diagonal  elements  of  H).  Thus  H  has  no  0  diagonal  entry,  resulting  in  an 
invertible  R  =  F~^HF.  The  only  issue  that  remains  now  is  of  the  conditioning  of  such  a  R.  The 
above  procedure  for  computing  R  results  in  a  relatively  well-conditioned  matrix,  with  the  condition 
number  of  R  ranging  from  24  for  iV  =  16  to  389  for  N  =  256.  Intermediate  values  of  N  result  in 
condition  numbers  between  24  and  389.  In  the  work  of  this  paper  we  use  this  practical,  and  thus 
invertible,  filtering  operator  given  in  (30)  for  all  calculations. 
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